---
title: "conjoint analysis"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
rm(list=ls())

require(tidyverse)
require(cjoint)
library(cregg)
library(scales)
library(patchwork)
library(extrafont)
require(memisc)

dat_legal <- read_csv("data/cj_data_legal.csv")
str(dat_legal)
dat_legal[,c(24:26)] <- data.frame(lapply(dat_legal[,c(24:26)],as.factor))

# H1: Islam
# Islamism
dat_legal <- dat_legal |> 
  mutate(Islamism = case_when(Q3_4 <= 2 ~ "High",
                              Q3_4 == 4 | Q3_4 == 5 ~ "Low",
                              Q3_4 == 3 | Q3_4 == 6 ~ "other"))
dat_legal$Islamism <- as.factor(dat_legal$Islamism)
#dat_legal$Islamism[dat_legal$Islamism == "other"] <- NA

dat_legal <- dat_legal |> 
  mutate(islamism = case_when(Q3_4 <= 2 ~ "High",
                              Q3_4 == 4 | Q3_4 == 5 | Q3_4 == 3 | Q3_4 == 6 ~ "Low"))
dat_legal$islamism <- as.factor(dat_legal$islamism)


# H2: Anti-U.S. by Shia parties supporter
dat_legal <- dat_legal |> 
  mutate(Shia_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Shia_party_supporter = ifelse(Shia_party_supporter >= 7, 1, 0))
dat_legal$Shia_party_supporter <- as.factor(dat_legal$Shia_party_supporter)


# H3: Liberalism 
# Liberalism as ideology
dat_legal <- dat_legal |>
  mutate(Liberalism = case_when(Q3_5 <= 2 ~ "High",
                                Q3_5 == 4 | Q3_5 == 5 ~ "Low",
                                Q3_5 == 3 | Q3_5 == 6 ~ "other"))
dat_legal$Liberalism <- as.factor(dat_legal$Liberalism)

# Liberalism by education
dat_legal <- dat_legal |> 
  mutate(Education = case_when(D3 <= 2 ~ "Low",
                               D3 == 3 | D3 == 4 ~ "Middle",
                               D3 >= 5 ~ "High"))
dat_legal$Education <- as.factor(dat_legal$Education)

# Liberalism by urban-ness
dat_legal <- dat_legal |> 
  mutate(Provincial_capital = case_when(D8 == 1 ~ "Urban",
                                        D8 == 2 ~ "Rural"))
dat_legal$Provincial_capital <- as.factor(dat_legal$Provincial_capital)

# Liberal_high_education
dat_legal <- dat_legal |>
  mutate(Liberal_high_education = ifelse(Liberalism == "High" & Education == "High", "High-educated_liberal", "Other"))
dat_legal$Liberal_high_education <- as.factor(dat_legal$Liberal_high_education)

# Liberal_urban
dat_legal <- dat_legal |>
  mutate(Liberal_capital = ifelse(Liberalism == "High" & Provincial_capital == "Urban", "Urban_liberal", "other"))
dat_legal$Liberal_capital <- as.factor(dat_legal$Liberal_capital)  


# H4: Conservative
# Conservative by tribalism
dat_legal <- dat_legal |>
  mutate(Tribalism = case_when(Q3_6 <= 2 ~ "Conservative",
                               Q3_6 == 4 | Q3_6 == 5 ~ "Non-conservative",
                               Q3_5 == 6 | Q3_6 == 6 ~ "other"))
dat_legal$Tribalism <- as.factor(dat_legal$Tribalism)

```

# conjoint setting
```{r}
f1 <- selected ~ amendment + reaction1 + reaction2
str(dat_legal)

dat_legal$amendment <- factor(dat_legal$amendment, 
                          levels = c("reform","against Freedom","Social Order"))
dat_legal$reaction1 <- factor(dat_legal$reaction1, 
                          levels = c("several reactions", "Western criticize", "Shia support"))
dat_legal$reaction2 <- factor(dat_legal$reaction2,
                            levels = c("commentator-one issue","sporkesperson-foreign interfere", 
                                       "professor-abusing human rights"))
```

# whole model
```{r, fig.height = 5, dpi = 300}
# MM
whole_mm <- cj(dat_legal, f1, 
               id = id,
               estimate = "mm", 
               h0 = 0.5)
whole_mm
plot(whole_mm, vline = 0.5, 
     size = 2)
write_html(whole_mm, file = "result/whole_mm.html")

# amce
whole_amce <- cj(data = dat_legal, f1,
                 id = ~ respondentIndex,
                 estimate = "amce")
whole_amce
plot(whole_amce, vline = 0,
     size = 2)
write_html(whole_amce, file = "result/whole_amce.html")

# plot estimation difference
whole_mm$Estimate <- "MM"
whole_amce$Estimate <- "AMCE"
plot(rbind(whole_mm, whole_amce), feature_headers = FALSE, size = 2)+
  ggplot2::facet_wrap(~ Estimate, ncol = 2L)
```

# H1: by Islamism
```{r, fig.height = 7, dpi = 300}
# amce
islam_amce <- cj(data = dat_legal,
                 selected ~ amendment + reaction1 + reaction2,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ Islamism)
islam_amce_diff <- cj(data = dat_legal,
                      selected ~ amendment + reaction1 + reaction2,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ Islamism)
plot(rbind(islam_amce, islam_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
dat_legal$Islamism <- relevel(dat_legal$Islamism, "Low")
islam_mm <- cj(dat_legal,
               selected ~ amendment + reaction1 + reaction2,
               id = id,
               estimate = "mm",
               by = ~ Islamism)
islam_mm_diff <- cj(data = dat_legal,
                    selected ~ amendment + reaction1 + reaction2,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ Islamism)
plot(rbind(islam_mm, islam_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# plot Islamism
```{r, fig.height = 5, dpi = 300}
plot(islam_amce, group = "Islamism")
plot(islam_amce, size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
plot(islam_mm, size = 2, vline = 0.5) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

plot(islam_amce_diff, size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
plot(islam_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~BY, ncol = 2L)
```

# OLS result of Islamism
```{r, fig.height = 5, dpi = 300}
#as.data.frame(islam_amce)
#as.data.frame(islam_mm_diff)
#as.data.frame(islam_amce_diff)

write_html(islam_amce, file = "result/islam_ame.html")
write_html(islam_mm, file = "result/islam_mm.html")
write_html(islam_amce_diff, file = "result/islam_ame_diff.html")
write_html(islam_mm_diff, file = "result/islam_mm_diff.html")

# anova
cj_anova(data = dat_legal, 
         selected ~ amendment + reaction1 + reaction2, by = ~ Islamism)
```

# H1: by Islamism2
```{r, fig.height = 7, dpi = 300}
# amce
islam_amce2 <- cj(data = dat_legal,
                 selected ~ amendment + reaction1 + reaction2,
                 id = ~ respondentIndex,
                 estimate = "amce",
                 by = ~ islamism)
islam_amce_diff2 <- cj(data = dat_legal,
                      selected ~ amendment + reaction1 + reaction2,
                      id = ~ respondentIndex,
                      estimate = "amce_diff",
                      by = ~ islamism)
plot(rbind(islam_amce2, islam_amce_diff2), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 2L)

# mm
dat_legal$islamism <- relevel(dat_legal$islamism, "Low")
islam_mm2 <- cj(dat_legal,
               selected ~ amendment + reaction1 + reaction2,
               id = id,
               estimate = "mm",
               by = ~ islamism)
islam_mm_diff2 <- cj(data = dat_legal,
                    selected ~ amendment + reaction1 + reaction2,
                    id = id,
                    estimate = "mm_diff",
                    by = ~ islamism)
plot(rbind(islam_mm2, islam_mm_diff2), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 2L)
```

# plot Islamism2
```{r, fig.height = 5, dpi = 300}
plot(islam_amce2, group = "islamism")
plot(islam_amce2, size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
plot(islam_mm2, size = 2, vline = 0.5) + 
  ggplot2::facet_wrap(~BY, ncol = 2L)

plot(islam_amce_diff2, size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 1L)
plot(islam_mm_diff2, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~BY, ncol = 2L)
```

# OLS result of Islamism
```{r, fig.height = 5, dpi = 300}
#as.data.frame(islam_amce)
#as.data.frame(islam_mm_diff)
#as.data.frame(islam_amce_diff)

write_html(islam_amce2, file = "result/islam_ame2.html")
write_html(islam_mm2, file = "result/islam_mm2.html")
write_html(islam_amce_diff2, file = "result/islam_ame_diff2.html")
write_html(islam_mm_diff2, file = "result/islam_mm_diff2.html")

# anova
cj_anova(data = dat_legal, 
         selected ~ amendment + reaction1 + reaction2, by = ~ islamism)
```


# H2: by Shia party supporter as indicator of anti-U.S.
```{r, fig.height = 5, dpi = 300}
# amce
shia_amce <- cj(data = dat_legal,
                selected ~ amendment + reaction1 + reaction2,
                id = ~ respondentIndex,
                estimate = "amce",
                by = ~ Shia_party_supporter)
shia_amce_diff <- cj(data = dat_legal,
                     selected ~ amendment + reaction1 + reaction2,
                     id = ~ respondentIndex,
                     estimate = "amce_diff",
                     by = ~ Shia_party_supporter)
plot(rbind(shia_amce, shia_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
shia_mm <- cj(data = dat_legal,
              selected ~ amendment + reaction1 + reaction2,
              id = id,
              estimate = "mm",
              by = ~ Shia_party_supporter)
shia_mm_diff <- cj(data = dat_legal,
                   selected ~ amendment + reaction1 + reaction2,
                   id = id,
                   estimate = "mm_diff",
                   by = ~ Shia_party_supporter)
plot(rbind(shia_mm, shia_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# anova
cj_anova(data = dat_legal, 
         selected ~ amendment + reaction1 + reaction2, by = ~ Shia_party_supporter)
```


# H3-1: by Liberalism
```{r, fig.height = 6, dpi = 300}
# amce
liberal_amce <- cj(data = dat_legal,
                   selected ~ amendment + reaction1 + reaction2,
                   id = ~ respondentIndex,
                   estimate = "amce",
                   by = ~ Liberalism)
liberal_amce_diff <- cj(data = dat_legal,
                        selected ~ amendment + reaction1 + reaction2,
                        id = ~ respondentIndex,
                        estimate = "amce_diff",
                        by = ~ Liberalism)
plot(rbind(liberal_amce, liberal_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
liberal_mm <- cj(dat_legal,
                 selected ~ amendment + reaction1 + reaction2,
                 id = id,
                 estimate = "mm",
                 by = ~ Liberalism)
liberal_mm_diff <- cj(data = dat_legal,
                      selected ~ amendment + reaction1 + reaction2,
                      id = id,
                      estimate = "mm_diff",
                      by = ~ Liberalism)
plot(rbind(liberal_mm, liberal_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# plot Liberalism
```{r, fig.height = 5, dpi = 300}
plot(liberal_amce, size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
plot(liberal_mm, size = 2, vline = 0.5) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

plot(liberal_amce_diff, size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
plot(liberal_mm_diff, size = 2, vline = 0) + 
  ggplot2::facet_wrap(~BY, ncol = 2L)
```

# OLS result of Liberalism
```{r, fig.height = 5, dpi = 300}
#as.data.frame(liberal_amce)
#as.data.frame(liberal_mm_diff)
#as.data.frame(liberal_amce_diff)

write_html(liberal_amce, file = "result/liberal_ame.html")
write_html(liberal_mm, file = "result/liberal_mm.html")
write_html(liberal_amce_diff, file = "result/liberal_ame_diff.html")
write_html(liberal_mm_diff, file = "result/liberal_mm_diff.html")

# anova
cj_anova(data = dat_legal, 
         selected ~ amendment + reaction1 + reaction2, by = ~ Liberalism)
```

# H3-2: by Liberalism and high education
```{r, fig.height = 5, dpi = 300}
# amce
liberal_hiedu_amce <- cj(data = dat_legal,
                        selected ~ amendment + reaction1 + reaction2,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ Liberal_high_education)
liberal_hiedu_amce_diff <- cj(data = dat_legal,
                              selected ~ amendment + reaction1 + reaction2,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ Liberal_high_education)
plot(rbind(liberal_hiedu_amce, liberal_hiedu_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
dat_legal$Liberal_high_education <- relevel(dat_legal$Liberal_high_education, "Other")

liberal_hiedu_mm <- cj(dat_legal,
                       selected ~ amendment + reaction1 + reaction2,
                       id = id,
                       estimate = "mm",
                       by = ~ Liberal_high_education)
liberal_hiedu_mm_diff <- cj(data = dat_legal,
                            selected ~ amendment + reaction1 + reaction2,
                            id = id,
                            estimate = "mm_diff",
                            by = ~ Liberal_high_education)
plot(rbind(liberal_hiedu_mm, liberal_hiedu_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)
```

# OLS result of Liberalism and high education
```{r, fig.height = 5, dpi = 300}
write_html(liberal_hiedu_amce, file = "result/liberal_hiedu_ame.html")
write_html(liberal_hiedu_mm, file = "result/liberal_hiedu_mm.html")
write_html(liberal_hiedu_amce_diff, file = "result/liberal_hiedu_ame_diff.html")
write_html(liberal_hiedu_mm_diff, file = "result/liberal_hiedu_mm_diff.html")

# anova
cj_anova(data = dat_legal, 
         selected ~ amendment + reaction1 + reaction2, by = ~ Liberal_high_education)
```

# H3-3: by Liberalism and urban resident
```{r, fig.height = 5, dpi = 300}
# amce
liberal_urban_amce <- cj(data = dat_legal,
                        selected ~ amendment + reaction1 + reaction2,
                        id = ~ respondentIndex,
                        estimate = "amce",
                        by = ~ Liberal_capital)
liberal_urban_amce_diff <- cj(data = dat_legal,
                              selected ~ amendment + reaction1 + reaction2,
                              id = ~ respondentIndex,
                              estimate = "amce_diff",
                              by = ~ Liberal_capital)
plot(rbind(liberal_urban_amce, liberal_urban_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
liberal_urban_mm <- cj(dat_legal,
                       selected ~ amendment + reaction1 + reaction2,
                       id = id,
                       estimate = "mm",
                       by = ~ Liberal_capital)
liberal_urban_mm_diff <- cj(data = dat_legal,
                            selected ~ amendment + reaction1 + reaction2,
                            id = id,
                            estimate = "mm_diff",
                            by = ~ Liberal_capital)
plot(rbind(liberal_urban_mm, liberal_urban_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# anova
cj_anova(data = dat_legal, 
         selected ~ amendment + reaction1 + reaction2, by = ~ Liberal_capital)
```


# H4: by Tribalism as indicator of conservative
```{r, fig.height = 6, dpi = 300}
# amce
tribalism_amce <- cj(data = dat_legal,
                     selected ~ amendment + reaction1 + reaction2,
                     id = ~ respondentIndex,
                     estimate = "amce",
                     by = ~ Tribalism)
tribalism_amce_diff <- cj(data = dat_legal,
                          selected ~ amendment + reaction1 + reaction2,
                          id = ~ respondentIndex,
                          estimate = "amce_diff",
                          by = ~ Tribalism)
plot(rbind(tribalism_amce, tribalism_amce_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# mm
tribalism_mm <- cj(dat_legal,
                   selected ~ amendment + reaction1 + reaction2,
                   id = id,
                   estimate = "mm",
                   by = ~ Tribalism)
tribalism_mm_diff <- cj(data = dat_legal,
                        selected ~ amendment + reaction1 + reaction2,
                        id = id,
                        estimate = "mm_diff",
                        by = ~ Tribalism)
plot(rbind(tribalism_mm, tribalism_mm_diff), size = 2) + 
  ggplot2::facet_wrap(~BY, ncol = 3L)

# anova
#cj_anova(data = dat_legal, 
#         selected ~ amendment + reaction1 + reaction2, by = ~ Tribalism)
```

# frequency and propotion
```{r}
plot(cj_freqs(dat_legal,
              f1,
              id = ~ respondentIndex))
```
